Method for modeling the transport of ions 
in hydrated cement systems 



field of the invention 



[01] The present invention relates to a computational tool that allows to evaluate the 
degradation of cement-based materials under various types of chemical attacks such as 
sulfates, chlorides, and plain water. It is based on the physical principles of ionic mass 
conservation and chemical equilibrium between a solution and different solid phases. 
Q The effect of the dissolution or the precipitation of solid phases on the transport 
|D coefficients is considered. The transport equations are discretized in a computer code 
m using the finite element method (FEM). 

O BACKGROUND OF THE INVENTION 

[02j^ The degradation of concrete structures under the effect of chemical attacks has become 
FU a major concern for civil engineers. The most common examples of chemical attacks 
O are the corrosion of reinforcement bars as a result of chloride exposure, external sulfate 
|a& attacks, carbonation and alcali-aggregate reactions. For the corrosion case, the cost of 

repair of the concrete structures exposed to this problem is estimated at 20 billion 

dollars in the US only. 



[03] All of these examples of premature concrete structure degradations have their origin in 
ionic exchanges between the material and its environment. The corrosion of 
reinforcement bars is caused mainly by the penetration of external chloride in the 
concrete porous network. The presence of chloride depassivates the steel, which 
initiates corrosion. The expansive corrosion product, when in sufficient amount, will 
cause the concrete to crack and eventually fall apart. External sulfate attack is based on 
the same principle. In that case, it is the penetration of sulfate ions in the concrete that 
is at the origin of the degradation. When the proper chemical conditions are met, the 
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sulfates ions react with the solid phases of the material to form gypsum and ettringite. If 
the amount of these two products is important, severe cracking can occur, again leading 
to a degradation of the structure. 

The reverse situation is also possible. Ions originally present in the concrete porous 
network as a result of cement hydration can be leached out of the material to the 
external environment. For example, a structure in contact with pure water will lose some 
of its hydroxyl ions, which causes a reduction of pH in the material. This drop of pH will 
lead to the dissolution of portlandite (calcium hydroxide) and the decalcification of C-S- 
H, the binding phase of the material. Consequently, the structure is bound to lose its 
mechanical resistance. 

A good knowledge of ionic transport mechanisms in cement-based materials and the 
implementation of the related physical laws in a computer model would be necessary to 
evaluate properly the service-life of concrete structures. Such a model could be used to 
plan reparation schedules for structures, based on their estimated remaining service- 
life. It could also help to choose the proper concrete mixture for a given structure, given 
the environment to which it is exposed. 

Despite this need, few significant developments related to the ionic transport in cement- 
based materials have been published in recent years. Most models used to predict the 
service-life of concrete structures are based on Fick's law, which simplifies the problem 
too much to yield reliable predictions. 

Fick's law accounts for the transport of ions as a result of diffusion. The transport of 
particles by diffusion is the result of their thermal agitation which causes random 
collisions that eventually disperses the particles from high concentration regions to 
weak concentrations regions. When the particles are electrically charged like ions, their 
charge influences the transport by diffusion through the electrical coupling between the 
ions and the chemical activity of the solution. As a consequence, the movement of one 
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ionic species is influenced by all the other species and the use of Fick's law is 
inappropriate. 

[08] The electrical coupling and the chemical activity effect emphasizes the multiionic aspect 
of ionic transport. Most of the time, the chloride ion is under scrutiny, as it is an 
important factor in the rebar corrosion phenomenon. However, few models consider the 
other ionic species involved in the transport process. They use Fick's law to model the 
transport of ions by diffusion. This approach is at the core of most ionic transport 
models in cement-based materials. 

[09J5 This is the case for example in the models published by Gospodinov et al. in "Diffusion 
of sulfate ions into cement stone regarding simultaneous chemical reactions and 
W resulting effects" published in "Cement and Concrete Research", vol. 29, p. 1591-1596 
O in 1999, by Hansen and Saouma, in "Numerical simulation of reinforced concrete 
^ deterioration - part 1: chloride diffusion", published in "ACI Materials Journal", vol. 96, 
Q no. 2, p.173-180, 1999, by Martin-Perez in "Service-life modeling of R.C. highway 
m structures exposed to chlorides", in a Ph.D. thesis, University of Toronto, 1998, by 
J Nagesh and Bhattacharjee in "Modeling of chloride diffusion in concrete and 
H determination of diffusion coefficients", published in ACI Materials Journal, vol. 95, no. 
2, p. 113-120, 1998, by Saetta et al. in "Analysis of chloride diffusion into partially 
saturated concrete", published in "ACI Materials" Journal, vol. 90, no. 5, p.441-451, 
1993, and Swaddiwudhipong et al. in, "Chloride ingress in partially and fully saturated 
concretes", published in "Concrete Science and Engineering", vol. 2, p. 17-31, 2000. 

[010] There are however an increasing number of papers where the electrical coupling is 
taken into account. This is the case for the work of Masi et al. in "Simulation of chloride 
penetration in cement-based materials", published in "Cement and Concrete Research", 
vol. 27, no. 10, p.1591-1601, 1997 and True et al. in, "Numerical simulation of multi- 
species diffusion", published in, "Materials and Structures", vol. 33, p.566-573, 2000. It 
should be emphasized that there are very few models accounting for the chemical 
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activity effects. Li and Page in, "Modeling of electrochemical chloride extraction from 
concrete: influence of ionic activity coefficients", published in "Computational Materials 
Science", vol. 9, p. 303-308, 1998 include it in their model, as well as the electrical 
coupling. However, their model is limited to ionic transport in saturated cement-based 
systems exposed to an electrical current. It is not relevant to predict the service-life of 
concrete structures. 

[011] Ions can also move under the effect of fluid displacement in the porous network of the 
material. This phenomenon is called advection. This fluid displacement occurs as a 
result of capillary forces arising in unsaturated materials. A structure becomes 
2 unsaturated following the various wetting/drying cycles to which it is exposed during its 
[2 service-life. In the models cited previously, only those of Martin-Perez, Nagesh and 
W Bhattacharjee, Saetta et aL, and Swaddiwudhipong et al. consider capillary forces in 
□ unsaturated cement-based materials as a transport mechanism. 

[012] Both diffusion and advection drag the ions through the porous network of concrete. The 
m ions may then undergo some chemical reactions with the hydrated cement paste. For 
;1 example, the penetration of sulfate (S0 4 2 ") ions in cement-based materials may lead to 
U the formation of ettringite and gypsum, while the penetration of chloride is at the basis of 
the formation of chloroaluminates. Studies performed on simple cement systems clearly 
showed that these chemical reactions are multiionic. For instance, in addition to SC^ 2 ", 
The formation of ettringite also involves Ca 2+ , OH" and AI(OH) 4 " ions. The last three ions 
are involved in the formation of chloroaluminates. Furthermore, the formation of 
ettringite, gypsum, and chloroaluminates is influenced by the presence of the alkali ions 
Na + and K + . 

[013] Some studies have indicated that surface binding phenomena (also called physical 
adsorption) can have a significant influence on ionic transport mechanisms. This 
appears to be particularly the case for chloride penetration problems. 
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[014] All of the previously cited models take into account interactions between ions in solution 
and the hydrated cement paste. They all use, without any exception, a simplified 
interaction model based on the concentration of only one ionic species. It is called the 
interaction isotherm, which consists in an experimental curve relating the amount of a 
given ion bound to the solid phase as a function of its concentration in solution. In most 
cases, the isotherm is determined for the chloride ions. This method is used in the 
model based on one ion as well as in the multiionic ones. Even if it allows to model the 
interactions involving chloride or sulfate, in most cases it does not allow to take into 
account other chemical reactions occurring simultaneously like the dissolution of 
portlandite or the decalcification of the C-S-H. The use of a simple chemical model thus 

2 limits the possibility of considering some reactions that are bound to have an important 
effect on the ionic transport in the material. Furthermore, it blends into one unique curve 

W the chemical reactions, where products are formed or dissolved, with the surface 

q interactions, where ions are adsorbed onto the surface. 

sssss 

[01 Therefore, since it is essential to be able to determine the behavior of hydrated cement 

ni systems and concrete structures, there is a need for a method of calculating such a 

£ service-life accurately. 

SUMMARY OF THE INVENTION 

[016] Accordingly, an object of the present invention is to provide a method and an apparatus 
for calculating a service-life for cement-based materials under chemical attacks. 

[017] Accordingly, another object of the present invention is to provide a calculation tool that 
allows to evaluate the variation over time of ionic profiles and solid phase contents in 
hydrated cement systems. It should take into account the following: the diffusion of 
multiple ionic species, the electrical coupling between ions, chemical activity effects, 
capillary suction, adsorption of ions and the dissolution/precipitation reactions. 
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[018] Accordingly, another object of the present invention is to provide a calculation tool that 
allows to evaluate the diffusion coefficient of ions in a cement-based material. 

[019] According to an aspect of the present invention, there is provided a method for 
determining an ion concentration in solution for each of at least two ions capable of 
undergoing transport in a cement-based material under a chemical attack and a solid 
phase profile of at least one component of the cement-based material. The cement- 
based material has a solid skeleton in which liquid-filled and/or vapor-filled pores are 
found. The porosity of the cement-based material is provided. The method comprises 
the steps of determining a first concentration for each ion and an electrical potential 
Jj profile using a transport algorithm, wherein the transport algorithm is a function of a 
52 diffusion of the ions, of an adsorption of the ions, of an electrical coupling between the 
y ions and of a chemical activity effect and wherein the ionic solution of the material is not 
in equilibrium with the various solid phases of the hydrated cement paste of the cement- 
+ based material; calculating a corrected concentration for each ion and a corrected solid 
O phase profile for the at least one of component using a chemical reactions algorithm, 
j|j wherein at least one of dissolution and precipitation reactions is accounted for in order 
*f to maintain an equilibrium between the ionic solution and the various solid phases of the 
U hydrated paste; calculating a changed transport properties profile to take into account 
an effect of the chemical reactions on the porosity of the material; and determining an 
ion concentration and a solid phase profile for at least one component using the 
changed transport properties profile and the corrected concentration and profile, 
wherein effects of the chemical reactions on the ionic transport are taken into account 
by correcting material transport properties and whereby the ion concentration for each 
ion and the solid phase profiles in the cement-based material can be used to evaluate a 
degradation of the cement-based material. 

[020] The transport algorithm can also be a function of the capillary suction if there is a 
relative humidity gradient. The water content is then further determined using said 
transport algorithm. 
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[021] According to another aspect of the present invention, there is provided a method for 
determining a diffusion coefficient for each of at least two ions capable of undergoing 
transport in a cement-based material, said cement-based material having a solid 
skeleton and pores, said pores being at least one of liquid-filled and vapor-filled, a 
porosity of said cement-based material being provided, the method comprising the 
steps of: determining a concentration for each said at least two ions and an electrical 
potential profile using a transport algorithm, wherein the transport algorithm is a function 
of a diffusion of said at least two ions, of an electrical coupling between said at least two 
ions and a chemical activity of each said at least two ions; and determining a diffusion 

*S coefficient for each of at least two ions using said concentration and said electrical 

J2 potential profile. 

[02|| Preferably, the method further comprises providing input data, said input data 
P comprising a set of tortuosity parameters and a measured electrical current for said 
o material, determining a diffusion coefficient for each of said ions by 1 . calculating a set 
Kj of electrical currents each electrical current corresponding to one of said tortuosity 
£ parameter in said set of tortuosity parameters using said concentration, said electrical 
M potential and each one of said tortuosity parameter in said set of tortuosity parameters; 
2. choosing an electrical current from said set with a value closest to said measured 
electrical current; 3. determining a matching tortuosity parameter corresponding to said 
chosen electrical current; 4. determining said diffusion coefficient for each of said ions 
using said matching tortuosity. 

[023] According to another aspect of the present invention, there is provided an apparatus for 
determining an ion concentration in solution of each of at least two ions capable of 
undergoing transport in a cement-based material under a chemical attack and a solid 
phase profile of at least one component of the cement-based material. The apparatus 
comprises a first determiner for determining a first concentration for each ion and an 
electrical potential profile using a transport algorithm; a first calculator for calculating a 
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corrected concentration for each ion and a corrected solid phase profile for the at least 
one of component using a chemical reactions algorithm; a second calculator for 
calculating a changed transport properties profile to take into account an effect of the 
chemical reactions on the porosity of the material; and a second determiner for 
determining an ion concentration and a solid phase profile for at least one component 
using the changed transport properties profile and the corrected concentration and 
profile. 

[024] According to still another aspect of the present invention, there is provided an apparatus 
for determining a diffusion coefficient for each of at least two ions capable of undergoing 
J transport in a cement-based material. The apparatus comprises a first determiner for 
m determining a concentration for each ion and an electrical potential profile using a 
% transport algorithm, and a second determiner for determining a diffusion coefficient for 
O each ion using the concentration and the electrical potential profile. 

[02- j Preferably, the apparatus further comprises an input data provider for providing input 
IVr data, the input data comprising a set of tortuosity parameters and a measured electrical 
O current for the material and a third determiner for determining a diffusion coefficient for 
^ each ion by 1. calculating a set of electrical currents each electrical current 
corresponding to one of the tortuosity parameter in the set of tortuosity parameters 
using the concentration, the electrical potential and each one of the tortuosity parameter 
in the set of tortuosity parameters; 2. choosing an electrical current from the set with a 
value closest to the measured electrical current; 3. determining a matching tortuosity 
parameter corresponding to the chosen electrical current; 4. determining the diffusion 
coefficient for each of the ions using the matching tortuosity. 

BRIEF DESCRIPTION OF THE DRAWINGS 

[026] These and other features, aspects and advantages of the present invention will become 
better understood with regard to the following description and accompanying drawings 
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wherein: 

[027] Figure 1 is a flow chart of an overview of a preferred algorithm according to the 
invention; 

[028] Figure 2a is an illustration of a representative elementary volume (REV); Figure 2b is an 
illustration of a plurality of REVs in which the pores are being filled with precipitates and 
solutes; 

[02&J Figure 3 illustrates a typical setup for the migration test; 

[03l] Figure 4 is a graph of the chloride concentration (mmol/L) versus the position (mm) in a 
;S sample problem. It gives the chloride concentration profiles after 10 hours when the 
O internal potential is either considered or not; 

[OM] Figure 5 is an illustration of the setup for the case considered for the numerical 
fU simulations; 

[ofe] Figure 6 is an illustration of the water content versus the relative humidity (%) for the 
adsorption isotherm according to Xi's model; 

[033] Figure 7 is an illustration of the concentration (g/kg) versus the position (cm) for the 
solid phase distribution after 20 years; 

[034] Figure 8 is an illustration of the concentration (mmol/L) versus the position in (cm) of a 
few ionic species and their distribution after 20 years; 

[035] Figure 9 is an illustration of the concentration (g/kg) versus the position (cm) for the 
comparison of the ettringite front after 10 years between the complete model and the 
cases where the electrical coupling and/or chemical activity effects are neglected; 
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[036] Figure 10 is an illustration of the concentration (g/kg) versus the position (cm) for the 
comparison of the ettringite and the portlandite profiles after 20 years for the saturated 
and unsaturated cases; 

[037] Figure 11 is an illustration of the concentration (g/kg) versus the position (cm) for the 
effect of the size of the time steps on the ettringite profile at 20 years; 

[038] Figure 12 is an illustration of the concentration (g/kg) versus the position (cm) for the 
,«! effect of the number of elements on the ettringite profile at 20 years; 

[03S] Figure 13 is an illustration of the experimental setup for the degradation of a cement 
^ paste exposed to Na2S0 4 solution. 

j sat 

[040] Figure 14 is an illustration of the current (mA) versus the time (hours) for the 

fi reproduction of the measured currents during the migration test with the second 

W embodiment of the invention; 

[oVl] Figure 15 is an illustration of the concentration (g/kg) versus the position (cm) for the 
solid phase distribution predicted by the model after 3 months of immersion in a 50 
mmol/L Na 2 S0 4 solution; 

[042] Figure 16 is an illustration of the calcium profiles after 3 months immersed in 50 mmol/L 
Na 2 S0 4 solution - The total calcium (g/kg) axis corresponds to the simulation and the 
count/second axis is associated to the microprobe data and both data sets are 
illustrated versus the position (mm); 

[043] Figure 17 is an illustration of the sulfur profile after 3 months immersed in 50 mmol/L 
Na2S0 4 solution - The total sulfur (g/kg) axis corresponds to the simulation and the 
count/second axis is associated to the microprobe measurements; 
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[044] Figure 18 is a flow chart of the main steps of the most preferred embodiment of the 
present invention; 

[045] Figure 19 is a flow chart of an overview of another preferred algorithm according to the 
invention; and 

[046] Figure 20 is a flow chart of the main steps of another preferred embodiment of the 
present invention. 

S DETAILED DESCRIPTION OF THE MORE PREFERRED EMBODIMENT 

[041^ An algorithm useful to the understanding of this invention was first described by Grove 
q and Wood in "Prediction and field verification of subsurface-water quality during artificial 
^ recharge", Lubbock, Texas, Groundwater, vol. 17, no. 3, p.250-257, 1979. It separates 
o the transport of ions from the chemical reactions. During a single time step, the 
equations for ionic transport are first solved, without considering any chemical reactions, 
jr At the end of this step, the concentration profiles are known. However, the 
U concentration of various points in the material may not be in equilibrium with the phases 
forming the solid skeleton of the material. Then starts the chemical step. The 
concentration profiles are corrected according to dissolution/precipitation reactions in 
order to maintain the equilibrium between the solution and the different solid phases. 
After the chemical step, the concentration profiles in solution as well as the solid 
concentration profiles are known. 

[048] The most preferred embodiment is summarized in Figure 1. The calculation starts from 
input data giving the initial state of the material, its transport properties and some 
numerical parameters: porosity, temperature, ionic pore solution concentrations, amount 
of each solid phase, diffusion coefficient for each ionic species, water diffusivity, initial 
water content, the number of time steps and their length, and the finite element mesh. 
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Boundary conditions are also provided to start the calculation. The calculation for the 
first time step begins. The invention performs a transport step, in which the ionic profiles 
are modified. These ionic profiles and the amount of the various solid phases are then 
corrected in the chemical step following dissolution and precipitation reactions in order 
to maintain the equilibrium between the pore solution and the hydrated paste. The 
invention then goes back to the transport step and begins the second time step, and so 
on. At the end of the calculation, the output of the invention is concentration and solid 
phase content profiles in the materials at various time steps. 

04&J It will be readily understood by one skilled in the art that in order to solve the transport 

3 equations, any numerical method could have been used. The preferred numerical 

S method for solving the transport equations is the finite elements method. However, the 

W finite differences and the finite volumes methods could be used without departing from 

O the present invention. 

06|1 The transport step 26 is summarized as follows. In the transport step, the movement of 
m ions is modeled by the extended Nernst-Planck equation with an advection term. This 
:t equation accounts for the electrical coupling between the ions as well as the chemical 
1=^ activity effects. The latter is calculated on the basis of a modified Davies equation that 
describes well the behavior of highly concentrated ionic solutions, such as those 
commonly found in cement-based materials. The advection term takes into account the 
movement of fluid arising when the material is exposed to gradient of relative humidity. 
A mass conservation equation is considered for each ionic species. This mass 
conservation equation includes a term to account for ionic adsorption. To estimate the 
electrical potential in the material caused by the electrical coupling between the ions, 
Poisson's equation is used. The last equation considered is Richards', which is used to 
calculate the water content distribution in the porous network. The system of equations 
is solved using the finite element method. 

[051] The chemical reaction model of dissolution and precipitation 27 is summarized as 
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follows. After a transport step, the solution in the pores is not in thermodynamic 
equilibrium with the various solid phases of the hydrated cement paste. To restore to 
the equilibrium state, the invention uses a chemical equilibrium code. It adjusts both the 
pore solution and the solid phase content. The equilibrium is calculated on the basis of 
the chemical equilibrium relationships of dissolution/precipitation reactions. It is 
assumed that the velocity of the ions is slow compared to the kinetics of the chemical 
reactions. This is called the local equilibrium hypothesis. Under this hypothesis, the 
chemical equilibrium relationships are valid. 

[052] During the course of the ionic transport process, the dissolution of some phases will 

J locally increase the porosity in the material, thus allowing the ions to move more rapidly 

JJ in a structure. On the contrary, the precipitation of some solid will reduce the section of 

y the pores and slow the ionic movement. The effects of the chemical reaction on the 

g ionic transport are taken into account by correcting the material transport properties as 

* a function of the porosity in step 28. This particular aspect of the invention is not present 

a in any of the prior art models. A check is made at step 29 to ensure that all time steps 

55; have occurred. The output 30 comprises ionic concentration profiles, solid phase 

•P profiles, water content profile and electrical potential profile. 

[053] A few assumptions had to be made prior to creating the model of the present invention. 
The model was designed under the hypothesis of a constant temperature. The solid 
phase is assumed non-deformable. External mechanical forces are thus not considered 
neither are cracks formed in the material. The calculations are performed assuming a 
constant hydration state. Finally, the hydrated products forming the solid phase are 
assumed to be uniformly distributed throughout the latter. The model is presented for 
the 1D case. Although these assumptions are used in the model of the preferred 
embodiment, it will be understood that modifications to the model to change these 
assumptions do not depart from the present invention. 

Transport model 
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[054] To model the transport of ions, which occur in the liquid (aqueous) phase, the equations 
are first written at the microscopic scale. They are then integrated over a 
Representative Elementary Volume (REV), as shown in Figure 2a, using the 
homogenization (averaging) technique, to yield the equations at the macroscopic scale. 
The REV comprises a solid skeleton 35 and pores containing solution 36 and/or a 
gaseous phase 37. An example of an arrangement of such pores is shown in Figure 2a. 
Details on the averaging technique can be found in Bear et al., "Introduction to Modeling 
of Transport Phenomena in Porous Media", Kluwer Academic Publishers (Netherlands), 
1991 and in Samson et al., "Describing ion diffusion mechanisms in cement-based 

J materials using the homogenization technique", Cement and Concrete Research, vol. 

j| 29, no. 8, p.1341-1345, 1999. 

[OfjJjJ The transport of ions may create precipitates 34 in the solution-filled pores 36 as is 
=P shown in Figure 2b. Figure 2b shows a plurality of REVs making up a portion of the 
P cement-based material. The precipitates 34, as the precipitation occurs, fill up the 
|[ solution-filled pores 36 and eliminate a portion of the solution-filled pores 36. This 
=F creates a change in the porosity of the material and affects its transport properties. The 
p solutes increase the porosity of the material by mixing in with the solution in the 

solution-filled pores 36 and increasing their size, see for example pore 36' which has 

been filled up by solutes. 

[056] The macroscopic equation for the transport of ionic species / is based on the extended 
Nernst-Planck equation with an advection term. Once integrated over the REV, it gives: 



afc c ;) ,d(fr,) a 

dt dt dt 



dx RT dx dx ~ 

V v ' v v ' ^ v ' advection 



y diffusion electrical coupling chemical activity 



= 0 (1) 



[057] where c, is the concentration of the species /' in solution, C/ s is the concentration in solid 
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phase, 0 S is the volumetric solid content of the material, 0 is the volumetric water 
content in the porous network, D, is the diffusion coefficient, z,* is the valence number of 
the species, F is the Faraday constant, R is the ideal gas constant, T is the temperature 
of the material, if/\s the electrical potential, y, is the chemical activity coefficient and V x is 
the average velocity of the fluid in the pore system under the action of capillary suction. 

[058] This equation is very different from classical diffusion model where only the term 
bearing the « diffusion » label in equation (1) usually appears. This term models the 
movement of the ions as a result of their thermal agitation. It is best known as Fick's 
law. 

[Ofjjtj] But ions in solution bear electrical charges. As ions have different drifting velocities, the 
W fastest ions tend to separate from the slower ones. However, since charges of opposite 
O signs mutually attract each other, the faster ions are slowed down and the slowest ones 
^ are accelerated, in order to bring the system near the electroneutral state. This creates 
J3 an electrical potential y/ in the material. This term is labeled « electrical coupling » in 
fy equation (1). To take this potential into account, Poisson's equation is added to the 
g model. It is given here after being averaged over the REV: 

ax\ ax J s J 

[060] where r is the tortuosity of the porous network, s is the dielectric permittivity of the 
solution and N is the total number of ionic species. The validity of this equation is based 
on the assumption that the electromagnetic signal travels much more rapidly than ions 
in solutions. This allows considering an equation from electrostatics in a transient 
problem, 

[061] The next term in equation (1) that differs from a classical diffusion model is related to 
the chemical activity of the ionic species. To calculate the chemical activity coefficients 
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yis, several models are available. However, well-known models such as Debye-Huckel 
or Davies are unable to cope with the specific case of cement-based materials, which 
bears a highly charged pore solution. A modification of Davies' relationship was found to 
yield good results: 



Azfyfl (0.2- 4.17 x\0~ 5 I)Azfl 
Iny, '■ — r + == '— (3) 

l + a ( .5V7 VlOOO 



[062] where / is the ionic strength of the solution: 



[063] In equation (4), A and B are temperature dependent parameters, given by: 



VlF 2 e 

A = — T» (5) 

%7i(eRTf' 2 



B = 



2F 2 



(6) 



sRT 

[064] where e 0 is the electrical charge of one electron, s=srs 0 is the permittivity of the medium, 
given by the dielectric constant times the permittivity of the vacuum. Finally, the 
parameter a, in equation (3) depends on the ionic species. Its value is 3x10" 10 for OH", 
3x10- 10 for Na + , 3.3X10" 10 for K + , IxlO' 10 for S0 4 2 ", 2x10" 10 for CI" and 1x10" 13 for Ca 2+ . 

[065] The movement of water in the pore network under the effect of capillary suction is 
modeled with Richard's equation. This equation is developed under the following 
hypothesis: isothermal conditions, isotropic material, non-deformable solid matrix, 
negligible gravitational effect, water movement slow enough to have equilibrium 
between the liquid and the gaseous phase, and a uniform pressure in both phases. The 



-16- 



15200-1 us JA/IC 



1 



equation is given by: 



dt dx 



D e — 
v fa J 



-0 (7) 



[066] where D 9 is the moisture diffusivity coefficient. It is the sum of the diffusivity coefficient 
of vapor and water. An expression is also needed for the average speed of the liquid 
phase, which appears in equation (1). It is given by: 



V.—dM (8) 

OX 



[Off] where D L is the water diffusivity coefficient. Equation (8) allows to transform equation (1 ) 

~ as: 



die A diec,) d( nr ,d Ci n D,z,F dy/ Qn d\n Yi d9^ 

V ■> 1 t j i Li fif) L _l fi — LJ — r — I— Jr. fin r — 4- H r 



■ + 

dt dt dt 



0D f — ^ + O^^c, + ODfii — + D L c, 

^ dx RT dx dx dx j 



(9) 



[0(38] The first term on the left-hand side of equation (9), involving c/\ describes the exchange 
between the solid and the aqueous phase, i.e. the ionic interactions. As shown in Rubin 
J., "Transport of reacting solutes in porous media: relation between mathematical nature 
of problem formulation and chemical nature of reactions", Water Resources Research, 
vol. 19, no. 5, p.1231-1252, 1983, this term can be used to model both the adsorption 
and the dissolution/precipitation reactions. For the model presented here, it will only be 
used to model adsorption. Adsorption is modeled through an adsorption isotherm that 
relates of to the concentration of that ionic species in solution c,-: Cj S =f(c } ). Depending on 
the material and on the ionic species in solution, the adsorption isotherm can have 
different forms. Here are some typical adsorption isotherm models: 

c- = Kc. linear isotherm (10) 
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< = Ac, 



Freundlich isotherm 



(ii) 



c- = 



Kc t 



linear isotherm 



(12) 



1 + Gc, 



where K, G and ah are constant to be determined experimentally. Once the isotherm is 
known, it is substituted in equation (9) knowing that: 



To get equation (14), it is assumed that the time variation of 0 S are small. 

The transport part of the model is complete. The unknowns are y/, #and N times c,, i.e. 
a concentration for each ionic species taken into account. There is accordingly N 
equations (14), one for each concentration, equation (2) to solve the electrical potential 
and equation (7) for the water content. The other parameters are either physical 
constants, like F, R, z- h e 0 , s, or parameters that have to be determined experimentally: 
T, e s , D u r, D 0} D Ll dof/dou 

To solve this system of nonlinear equations, the finite element method is used. Details 
are given in the following paragraphs. To ease the writing, the numerical model will be 
shown for a two ion case in 1D. But the model can easily be extended to take into 
account more species. The numerical examples discussed below involve six ions. 



dc; dc; dc, 



(13) 



dt dC; dt 



which gives : 




= 0 (14) 
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V 



[074] The finite element method requires one to write the equations under their weighted 
residual form. The weighted residual form, upon integration by parts, is given as: 



dc x 



0 
0 



dc 2 

0 
0 



c, 0 



c 2 0 



1 0 
0 0 



9 



dx 



dx 



9D X 

0 
0 
0 



0 

9D 2 
0 
0 



D L c 2 



9c x 



RT 

RT 2 
0 

6x 



C 2,x 



dx 



+ 



6D t 



1W 



dx 

0 

0 
0 

0 
0 
0 

--2,9 

s 



0 




0 


0 




0 


0 < 


dx 








0 




0 


0 


0 




0 


oj 


0 


0 


0~ 


V 


0 


0 


0 




0 


0 


0 


< ) 
9 


F 


0 


0 


z 2 9 

s 





9 



dx 



>dx-0 



(15) 



[075] where {8W) is the vector of the weighting function, defined as: 

= Sc x 89 Stff) (16) 



[076] The boundary terms are omitted since only Dirichlet and natural boundary conditions 
will be considered in the simulations. 
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[077] The weak form is discretized using the well-known Galerkin method with a standard 
linear two-node element (see Reddy J.N., Gartling D.K., "The finite element method in 
heat transfer and fluid dynamics", CRC Press (USA), 1994). for a complete text on the 
method). The approximation of the solution on each element is written as : 



6 
[N] = 



(17) 





0 


0 


0 


N 2 


0 


0 


0 


0 




0 


0 


0 


N 2 


0 


0 


0 


0 




0 


0 


0 


N 2 


0 


_ 0 


0 


0 




0 


0 


0 






C 2l 


0. 




C \2 


c 22 


°2 


¥2) 



(18) 



(19) 



where Ni and N 2 are the shape functions. The subscripts / and j for the concentrations 
in the vector (U n ) (equation 19) designate the species / at the node j of one element. 
The elementary matrices are thus expressed as : 



(20) 



M= lilNf[clN]) dx 



(21) 



with 
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N u 


0 


0 


0 


N 2, 


0 


0 


0 


0 


N u 


0 


0 


0 




0 


0 


0 


0 




0 


0 


0 




0 


0 


0 


0 




0 


0 


0 





(22) 



6Dj 


0 


D L c x 


0 


6D 2 


D L c 2 


0 


0 


D e 


0 


0 


0 



6bj 



RT 

o 

0r 



(23) 



[D 2 } = 



6D, 



0 
0 
0 



0D. 



0 

d(lny 2 ) 



dx 



0 
0 



0 0 

0 0 

0 0 

0 0 



(24) 



[A] = 



0 
0 
0 



F_ 

£ 



0 

0 



z x 6 



F_ 
e 



0 0 
0 0 
0 0 



z 2 0 0 0 



(25) 



[c]= 



0 + 0.** 



5 9cj 



0 
0 



0 + 0*'* 



0 
0 



c, 0 



c 2 0 

1 0 
0 0 



(26) 



[079] The assembly of elementary matrices [K®] and [M 6 ] leads to the following system of 
equations: 
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[Mfl}+[K]{U} = 0 (27) 



The time discretization is performed using the implicit Euler scheme: 
[M,]p^} + [*,]{C/,} = 0 (28) 

where At is the time step. The subscript t stands for the actual time step and t-At the 
previous one. Defining the matrices 

[fl-fMj+Afc] (29) 

[fHmJU^J (30) 
the system of equations can be written as: 

(31) 

where {Ut} is the vector of unknowns. This nonlinear system of equations is solved at 
each time step with the modified iterative Newton-Raphson method. Numerical 
simulations have shown that the convergence rate with a tangent matrix calculated 
without the nonlinear terms arising from the coupling between the various variables in 
the model is almost as fast as the one with a complete tangent matrix. However, the 
calculation time is reduced since less terms need to be calculated. Furthermore, its 
wider radius of convergence makes this algorithm even more interesting. The 
elementary tangent matrix is thus given as: 



[^]=[m 6 ']+A^] 



(32) 
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[084] Following Newton's method, the solution {U$ is obtained by successive iterations k 
according to : 

fc-']{A^}=-([^ t - 1 k- , }-{^- 1 }) (33) 

[085] This system of equations is solved for {Ai/} with Gauss elimination method and the 
solution is updated at each iteration : 

□ {U!}=\ur}+M (34) 

[Oil] Several iterations are performed over equations (33) and (34) until the increment {Alf} 
=fi becomes sufficiently small. When convergence is reached, the solution {U t } at this time 
J step is known, i.e. the variables c,, ^and 0are known. 

[081] The elementary matrices in equations (20) and (21) are evaluated through Gauss 
;| numerical integration method. Accordingly, the variables appearing in the matrices [Di], 
;=f [D 2 ], [D 3 ], and [C] are calculated at the integration points. To calculate the terms 
8{\ny,)/dx in (24), the ionic strength (equation 4) is calculated at each node. Then the 
Iny/S are calculated for each species and at each node with equation (3). From there, the 
value of 5(lny,)/5x can be evaluated at every integration point. 

Chemical reaction model 

[088] The chemical reactions considered in the invention are of the dissolution/precipitation 
type. The solids considered are listed in table 1. However, the algorithm presented is 
very general and could be applied to other solid phases as well. 
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[089] Table 1 : Solid phases in hydrated cement paste 



Name 


Chemical formula 


Expression for equilibrium 


Value* 


Portland ite 


Ca(OH) 2 


K sp ={Ca}{OHr 


5.2 


C-S-H 


1.65CaO.Si0 2 .(2.45)H 2 0 § 


K sp ={Ca}{OH} 211 


5.6* 


Ettringite 


3CaO.AI 2 03.3CaS04.32H 2 0 


K S p={Ca} 6 {OH} 4 {S0 4 } 3 {AI(OH)4} 2 


44.4 


Friedel's salt 


3CaO.AI 2 O 3 .CaCI 2 .10H 2 O 


K S p={Ca} 4 {OH} 4 {CI} 2 {AI(OH) 4 } 2 


29.1 


Hydrogarnet 


3CaO.AI 2 0 3 .6H 2 0 


K sp ={Ca} 3 {OH} 4 {AI(OH) 4 } 2 


23.2 


Gypsum 


CaS0 4 .2H 2 0 


K sp ={Ca}{S0 4 } 


4.6 


Mirabilite 


Na 2 SO 4 .10H 2 O 


K sp ={Na} 2 {S0 4 } 


1.2 



[098 Where {...} indicates chemical activity, § indicates that the C-S-H is considered having a 

S CIS ratio of 1 .65, 11 indicates that the C-S-H decalcification is modeled as the portlandite 

? S dissolution with a lower solubility and # indicates the value of the equilibrium constant (- 

O log Ks P ). 

[093] The equilibrium between the solid phases and the solution is checked at each node of 

ffi the finite element mesh. It is calculated as follow. Consider for example the case of 

S portlandite (see Table 1). Its equilibrium constant is expressed as: 

K sp ={Ca}{OHf (35) 

[092] where the curly brackets {...} indicate chemical activity. Knowing that the activity is equal 
to yfii and that in portlandite, for each Ca 2+ there is two OH", equation (35) is 
transformed as: 

+ Ac Cff fc„+2Ac c J (36) 

[093] where the subscript o indicates the concentration value before equilibrium. This 
equation must be solved for Ac Ca , which gives the correction needed to reach 
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[089] Table 1 : Solid phases in hydrated cement paste 



Name 


Chemical formula 


Expression for equilibrium 


Value* 


Portlandite 


Ca(OH) 2 


K sp ={Ca}{OH}' 


5.2 


C-S-H 


1.65CaO.Si0 2 .(2.45)H 2 0 § 


K sp ={Ca}{OH} 211 


5.6 11 


Ettringite 


3CaO.AI 2 0 3 .3CaS0 4 .32H 2 0 


K sp ={Ca} 6 {OH} 4 {S0 4 } 3 {AI(OH)4} 2 


44.4 


Friedel's salt 


3CaO.AI 2 0 3 .CaCI 2 .1 0H 2 O 


K sp ={Ca} 4 {OH} 4 {CI} 2 {AI(OH) 4 } 2 


29.1 


Hydrogarnet 


3CaO.AI 2 0 3 .6H 2 0 


K sp ={Ca} 3 {OH} 4 {AI(OH) 4 } 2 


23.2 


Gypsum 


CaS0 4 .2H 2 0 


K sp ={Ca}{S0 4 } 


4.6 


Mirabilite 


Na 2 SO 4 .10H 2 O 


K sp ={Na} 2 {S0 4 } 


1.2 



[09J] Where {...} indicates chemical activity, § indicates that the C-S-H is considered having a 

;= C/S ratio of 1.65, 11 indicates that the C-S-H decalcification is modeled as the portlandite 

W dissolution with a lower solubility and * indicates the value of the equilibrium constant (- 

O log Ksp)- 

[0S|] The equilibrium between the solid phases and the solution is checked at each node of 

ill the finite element mesh. It is calculated as follow. Consider for example the case of 

:jr portlandite (see Table 1). Its equilibrium constant is expressed as: 

K sp = {Ca}{OHf (35) 

[092] where the curly brackets {...} indicate chemical activity. Knowing that the activity is equal 
to rid and that in portlandite, for each Ca 2+ there is two OH", equation (35) is 
transformed as: 

K sp = YcaYln [tea + AC Ca \c° OH + l^Ca )' (36) 

[093] where the subscript o indicates the concentration value before equilibrium. This 
equation must be solved for Acc a , which gives the correction needed to reach 
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equilibrium. This procedure must be applied to all solid phases. At the end of this 
calculation step, all solid phases are in equilibrium with all ions in solution. The 
algorithm is detailed here : 

[094] 1 . Calculate the activity product ycayonCcaOoH ■ 

[095] 2. If the activity product is under the equilibrium constant and there is still some 
portlandite left, then there will be dissolution. If the activity product is over the 
equilibrium constant, there will be precipitation. In both cases, equation (36) is solved 
according to these calculation steps: 

[OQI] a) Solve equation (36) assuming constant chemical activity coefficients with the 
W Newton-Raphson method. This gives the following recursive formula: 

T k _ M r fa r 2 oH (c° Ca + Ac*; 1 fc„ + 2 ac£ ) 2 - k sp 

[Op] As stated earlier, the subscript o indicates the concentration value before the 

O equilibrium calculation. 

[098] b) Calculate new chemical activity coefficients with the concentrations (c C a 0+ Ac C a k ) 
and (c 0 H 0+ 2Acca k )- 

[099] c) Go back to step (a) with the new chemical activity coefficients. 

[0100] d) Repeat steps (a) to (c) until the activity product yc a yoH 2 cc a coH is near the 
equilibrium value K sp , according to a given threshold value. When equilibrium is reached 
for that particular phase, the concentrations are corrected with the last value of Ac Ca k 
found. 



-25- 



15200-1 us JA/IC 



[0101] 



3. Repeat steps (1) and (2) for the other solid phases (see table 1). 



[0102] 4. After each solid phase has been considered, verify if the corrected solution is in 
equilibrium with all solid phases. If this is not the case, repeat steps (1) to (3) until the 
equilibrium state is reached (again within a given threshold value). 

[0103] After this calculation, the concentrations in solution have been adjusted in order to 
respect the equilibrium expressions given in Table 1, indicating that dissolution and/or 
precipitation reactions occurred. The solid phases must be corrected accordingly at 
each node of the finite element mesh. To do so, the variations in concentration for each 
3 ions must be associated with the proper solid. For example, the total variation of 
S calcium, before and after equilibrium was made, can come from the dissolution of 
W portlandite, the dissolution of hydrogarnet and the precipitation of ettringite, since all 
5 these solids contain calcium. Assume that at a given node, there are variations of 
f concentration in Ca 2+ , S0 4 2 " and AI(OH) 4 ", and that portlandite and hydrogarnet are 
O present at this location. One can thus write : 



\Tot _ KEttringile /^ON 
A S04 _ a S04 v JO -' 



pTol _ ^Ettringite ^Hydrogarnet (39) 

£ot _ ^Ettringite ^Hydrogarnet _j_ ^Portlandite (40) 

[0104] where the A Tot quantities correspond to the variation in concentration after and before 
the equilibrium was made. The other zls correspond the fraction of the variations 
associated with each solid involved. Equations (38) to (40) are solved, knowing that 
AAI £ttnngite _ 2/3 A S o4 Ettrin9ite since in ettringite, for 2 moles of alumina, there is 3 moles of 
sulfate, and that A Ca Ettr,n9ite = 2A S o4 Btrin9ite and A Ca Hydr ° 9arnet = 3/2 A AI Hydf09amei \ for similar 
reasons. Once A S04 Emngite , A A , Hydr09arnet , and A Ca Portlandite are know, they are translated 
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into the proper amount of solid. 

[0105] The study of cement systems exposed to sulfate ions revealed that hydrogarnet and 
gypsum are two phases that do not coexist. Gypsum can only form after all hydrogarnet 
has been dissolved. These observations are taken into account in step 3, to avoid 
performing equilibrium on both hydrogarnet and gypsum. 

[0106] The C-S-H is treated in a particular way. When in contact with a low pH solution, the C- 
S-H phase loses its calcium, a phenomenon called decalcification. It leaves in place a 
silica gel. To model this chemical reaction, the C-S-H will be considered as portlandite, 
S but with a lower solubility constant (see Table 1). Its decalcification will thus release 
5™ Ca 2+ and OH" ions. 

[0|§7] The chemical reactions, besides binding or releasing ions while solid phases are 
* precipitated or dissolved, will have an effect on the transport properties by affecting the 
Q porosity of the material. If, for example, gypsum is formed at some place, it will reduce 
the porosity at that location, thus reducing the section across which ions are able to 
;| diffuse. This will reflect on their diffusion coefficient, according to an equation presented 
M by Garboczi and Bentz in "Computer simulation of the diffusivity of cement-based 
materials", Journal of Material Science, vol. 27, p.2083-2092, 1992: 



[0108] where <j> ca p is the capillary porosity of the paste, Df is the diffusion coefficient of ionic 
species / in free water (as opposed to a porous material) and H is the Heaviside function 
such that H(x)=1 for x>0 and H(x)=0 for x<0. The initial capillary porosity can be 
calculated with the following relationship: 




(41) 
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nm, (w/c)- 0.36a 

9cm = ■ 



(w/c) + 032 



(42) 



[0109] where w/c is the water/cement ratio of the cement paste and a is the degree of 
hydration of the paste (0 <a<1). Upon chemical modification to the paste, the capillary 
porosity is calculated as: 



(43) 



where V s is the volume of a given solid phase, per unit volume of cement paste, and M 
is the total number of solid phases. Based on this model, the correction factor G that 
multiplies the diffusion coefficients D-, of each ionic species is calculated as follow: 



Dt 



G = 



Modified paste 



D, 



Initial paste 



(44) 



[0111] Again, the case of C-S-H is treated separately. The silica gel remaining after the 
decalcification process has no structural strength. When at a given location the 
decalcification is completed, the damage factor G is set to ten. 



[0112] Figure 18 is a flow chart of the main steps of the most preferred embodiment of the 
present invention. The preferred embodiment is a method for determining an ion 
concentration for each of at least two ions capable of undergoing transport in a cement- 
based material under a chemical attack and a solid phase profile for at east one 
component of the cement-based material. The method comprises the steps of 
determining 70 a first concentration for each ion and an electrical potential profile using 
a transport algorithm. The water content profile can also be determined by the transport 
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algorithm. The next step is calculating 71 a corrected concentration for each ion and 
solid phase profiles using a chemical reactions algorithm, wherein dissolution and/or 
precipitation reactions are accounted for. The next step is calculating 72 a changed 
transport properties profile to take into account an effect of the chemical reaction 
equilibrium on the porosity of the material. Finally, the last step is determining 73 an ion 
concentration and solid phase profiles using the changed transport properties profile 
and the corrected concentration and profile. 

DESCRIPTION OF ANOTHER PREFERRED EMBODIMENT 

[0ffl3] The main parameter that characterizes ionic transport in porous materials is the 
^ diffusion coefficient D h This parameter can be determined with diffusion tests, but they 
yj are time consuming, their duration extending over several months. Since two decades, 
q the migration test is more commonly used. The experimental setup is similar to the 
* diffusion test (see Figure 3). It consists in a cement-based material disk 40 located 
o between two cells containing different solutions. Both cells contain a solution of NaOH 
Ji ! in order to avoid the degradation of the material. The upstream cell 38 also contains one 
:£ of these salt: NaCI, KCI, Na2S04. Prior to the test, the disk is immersed in a NaOH 
lI solution under vacuum for at least 24 hours to saturate the pores. Watertight joints 
prevent any leakage between the sample and the cells. As a result of diffusion, the ions 
of this salt will move through the sample to the downstream cell 41. The measurement 
of the ionic flux through the sample is used to evaluate the diffusion coefficients of the 
ionic species. In the case of the migration test, the ions are accelerated with an 
electrical field 43, which reduces the duration of the test considerably. 

[0114] The method for modeling presented herewith can be used to model the migration test. 
For the specific case of this test, some simplifications to the mathematical equations 
can be made. The migration tests are performed in saturated conditions, and no 
pressure gradient is applied on the liquid phase. In that case, the advection term in 
equation (14) can be dropped, as well as Richard ! s equation (equation 7). Furthermore, 
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the saturated conditions allow to write: 

9 = <f> (45) 
0,=1-* (46) 

[0115] where <j> is the porosity of the material. Finally, the short duration of the test (about 120 
hours), combined with the high velocity of the ions under the applied voltage, allow to 
neglect the chemical reactions between the ions and the hydrated paste of the material. 
Accordingly, the mass conservation equation (14) reduces to the extended Nernst- 
O Planck equation: 

w a A _ l( D , ^ + ^£ c , + Dfit ^0 = 0 (47) 
q dt dt\ ' dx RT dx ' ' 8x ) 

[01=46] Poisson's equation (2) is simplified to: 

ill! i 

a ^ + 7l§H = ° <48) 

[01 1 7] The flux of ion in the material is given as: 

ox RT ox ox 
[0118] In equation (48) and (49), the diffusion coefficient D, is defined as: 

D,=rD? (50) 

[0119] where r is the tortuosity of the material and Df is the diffusion coefficient of the species / 
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in free water, of values which can be found in physics handbooks. For example, the 
values of D, M for the most common species found in cement-based materials are: 
5.273x1CT 9 m 2 /s for OH", 1.334x1fJ 9 m 2 /s for Na + , 1.957x10" 9 m 2 /s for K + , 1.065x10" 9 
m 2 /s for S0 4 2 ", 0.792x10" 9 m 2 /s for Ca 2+ , and 2.032x10" 9 m 2 /s for CI". These values are 
constant and represent the diffusion coefficients in very dilute conditions. The tortuosity 
appears in the model as a result of the averaging procedure. It accounts for the intricate 
path the ions travel through in a porous material. Equation (50) has very important 
implications. It shows, for instance, that if the diffusion coefficient D, of one ionic species 
is known, t is also known and accordingly, the diffusion coefficient of the other ionic 

n species are known. It also shows that as long as the tortuosity of the material remains 

;5 unchanged, the D,s remain constant. 

[OigO] Equation (47) is commonly used to model the migration test. But no reliable numerical 
O solutions were available to solve it until recently. Accordingly, several simplifying 
7 assumptions have been used to solve this equation in order to analyze the results of 
migration tests. Two major trends are found in the literature: steady-state and non- 
Hj steady-state analysis. 

[Oftl] The assumptions underlying the steady-state analysis are given as follows. The steady- 
state is reached when the measured amount of chloride ions in the downstream cell 
varies linearly through time. This indicates a constant flux, which is the basic definition 
of the steady-state. The steady-state needs a very long time to be reached in classical 
diffusion test, but when the ions are accelerated with an external voltage, it is reached in 
a few days. During steady-state, all terms related to time in equation (47) are set to 
zero. This is equivalent to solving the flux equation (equation 49) with J; as a constant. 

[0122] To further simplify the analysis, many authors have also assumed that the external 
voltage is strong enough to neglect the internal potential arising from the difference in 
the drift velocity of the various ionic species. This hypothesis allows to write: 
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d JL = ^Im. = _£ = constant (5 1) 

dx L 

[0123] where Ay/ ex t is the applied voltage difference across the sample, L is the thickness of the 
specimen, and E ex t is the corresponding constant electric field. It is also assumed that 
the applied voltage is strong enough to neglect diffusion. Finally, the chemical activity is 
dropped from the analysis. In that case, equation (49) reduces to: 

j. = ^£j^L Cl ^m (52) 

RT L 

[0if4] where J-, is a constant since steady-state is reached. 

[0f|5] Another way to analyze the migration test is through a non-steady-state (or transient) 
=P analysis. In this case, one of the main assumption is that the chemical reactions can be 
□ neglected. Tests performed in transient state are rather short and, combined with the 
S high velocity of the ions under the applied field, the chemical reactions allegedly do not 
+ have time to have a significant effect on the ionic transport. Neglecting the chemical 
!l reactions also implies that there is no change to the microstructure of the paste during 
the duration of the test, which allows to consider the porosity and tortuosity as 
constants. As for the steady-state analysis, the external potential is assumed to be 
much stronger than its internal counterpart. The effect of the chemical activity is also 
neglected. With all these hypothesis, equation (47) reduces to: 

^l( Di d JL + W CiE ) = 0 (53) 
3t dt\ 1 dx RT ") 

[0126] and can be solved analytically. For semi-infinite media, the analytical solution of 
equation (53) with a constant D-, is: 
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2 



erfc 



K 24D t t 



2RT 



+ e 



1 RT Urfc 



V 



2JD,t 



2RT 



(54) 



[0127] where c 0 is the boundary condition at x=0. To obtain the diffusion coefficients with this 
model, a D, is fitted to chloride profiles measured on the samples undergoing the 
migration test. 

[0128] Both methods (steady-state and transient), are based on the hypothesis that the 
electrical coupling between the ions can be neglected when an external potential is 

0 applied. To verify the validity of this hypothesis, a sample problem is considered. The 
!| problem will be solved with and without the constant field hypothesis to see if the latter 
j 1 ! can be used in the analysis of the migration test. A cement-based material with a 
5 tortuosity of z=1/100 and a porosity of 30% is assumed. According to equation (50), the 
3 diffusion coefficient of each ionic species is: 5.273x1 0" 11 m 2 /s for OH", 1 .334x1 0" 11 for 
I a Na + , 1.957x10" 11 for K + , 1.065X10' 11 for S0 4 2 -, 0.792x10- 11 for Ca 2+ , and 2.032x10" 11 for 
J CI". The sample has a 25 mm thickness and is subjected to a 10V applied potential (this 

1 corresponds to a 400 Volt/m electrical field). The boundary conditions correspond to the 

2 typical setup as shown in Figure 3, i.e. OH=300, Na + =800, and Cl=500 mmol/L at x=0, 
and OH=300 and Na + =300 mmol/L at x=25 mm. All of the other concentrations are set 
to zero on both ends of the material. Finally, the initial concentrations in the sample are: 
OH =354.0, Na + =250.0, K + =120.0, SO 4 2 "=10.0, Ca 2+ =2.0 and Cl"=0.0 mmol/L. The initial 
potential is zero everywhere in the material. A temperature of 25°C is assumed. The 
chloride concentration profiles after 10 hours are shown in Figure 4. The huge 
difference between the solution of the extended Nernst-Planck model and the model 
with the constant electrical field shows that the fact of neglecting the internal potential 
when modeling a migration test can lead to very misleading results with regard to the 
chloride profiles, and accordingly to the ionic flux evaluation (see equation 49). 

[0129] The latter results lead to the development of a new method to analyze the migration 
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test. The new method proposed to evaluate the diffusion coefficients is based on non 
steady-state current measurements. The boundary conditions correspond to the 
concentrations in both cells as well as the imposed potential difference over the sample. 
The initial conditions are given by the concentrations of the ions in pore solution prior to 
the test. With these data, the equations are solved with different tortuosity values 
according to the numerical method presented in the preceding section. The numerical 
current l c num is calculated at the measurement times according to : 

/« = S FY j z,J i (55) 

[0-jlo] where S is the exposed surface of the sample and J, is the ionic flux given by equation 
W (49). For each tortuosity value, the error between the model and the measurements is 
O calculated as: 



jt(Q S -KT) 2 ( 56 ) 

[Oftl] where M is the total number of measurement and l c mes and I c nm are the measured and 
numerical currents. The tortuosity value leading to the smallest error with the 
measurements gives the proper diffusion coefficients for each ionic species in the 
material considered. 

[0132] This other preferred embodiment is summarized in Figure 19. The calculation starts 
from input data giving the material properties, numerical parameters and boundary 
conditions 75. The material properties can be the porosity, the temperature, ionic pore 
solution. The boundary conditions can be the concentrations and the electrical potential. 
The numerical parameters can be the number of time steps and their length, the finite 
element mesh and the set of tortuosity. 
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error = 



[0133] The invention first performs a transport step 76, in which the concentration profile for 
each ion and the electrical potential profile are determined using the diffusion 
coefficients corresponding to the first tortuosity. 

[0134] It will be readily understood by one skilled in the art that in order to solve the transport 
equations, any numerical method could have been used. The preferred numerical 
method for solving the transport equations is the finite elements method. However, the 
finite differences and the finite volumes methods could be used without departing from 
the present invention. 

[01§5] Then, a step of calculation of the current 77 follows. The electrical current is determined 
for each tortuosity. 

[0156] The transport step 76 and the current calculation step 77 are repeated for each 
* tortuosity parameter until a numerical current matching a measured current is found 
O which will give the diffusion coefficients for the material. 

[0j§7] Figure 20 is a flow chart of the main steps of this preferred embodiment of the present 
C invention. The preferred embodiment is a method for determining an ion diffusion 
coefficient for each of at least two ions capable of undergoing transport in a cement- 
based material. The method comprising the steps of: determining a concentration 80 for 
each the ions and an electrical potential profile using a transport algorithm and 
determining a diffusion coefficient 81 for each ion using the concentration and the 
electrical potential profile. 

[0138] Preferably, the method further comprises calculating 82 a set of electrical currents, each 
electrical current corresponding to one tortuosity parameter using the concentration, the 
electrical potential and the tortuosity parameter; choosing 83 an electrical current from 
the set with a value closest to the measured electrical current and determining a 
matching tortuosity parameter corresponding to the chosen electrical current and 
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determining 84 the diffusion coefficient for each ion using the matching tortuosity. 
FIRST EXAMPLE 
Description of the problem 

[0139] The first example concerns the prediction of the behavior of a concrete slab undergoing 
external sulfate attack. This type of attack is characterized by the penetration of sulfate 
(S0 4 2 ") in the material, which precipitates into ettringite and gypsum. The formation of 
these two phases leads to microcracking of the cement paste. The invention allows to 
S locate the ettringite and gypsum through time, giving valuable information on the 
1 degradation state of a structure. Parallel to this, the leaching of OH" and Ca 2+ out of the 
W material leads to the dissolution of portlandite and the decalcification of C-S-H, which 
3 weaken the structure. Again, the invention allows to track the dissolution fronts of both 
of these phases. 

[0M0] The case considered for the numerical simulations is shown in Figure 5. It consists of a 
£ 15 cm thick concrete slab S1 lying on a soil S2 bearing a high concentration of sodium 
S sulfate. The concrete S1 has a water/cement ratio of 0.50 and is made with an ASTM 
type V cement. The soil S2 on which the slab lies is at a relative humidity of 85%. The 
external environment 50, in contact with the upper part of the slab, is at a relative 
humidity of 70%. This gradient in relative humidity will create a flow of water directed 
upward that will contribute, along with diffusion, to move the S0 4 2 " ions into the material. 
Besides S0 4 2 ", five other ions are considered in the simulations: OH", Na + , K + , Ca 2+ , and 
AI(OH) 4 ". 

[0141] The boundary conditions are as follows. On the lower part of the slab S2 (x=0), a 
concentration of 5000 ppm of S0 4 2 " is imposed, which corresponds roughly to 50 
mmol/L. Accordingly, a concentration of 100 mmol/L of Na + is set at x=0. The 
concentration of all other species at x=0 is zero. On the upper part of the slab S2 
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(x=0.15), there is no transfer of ions from the material to the environment. This 
boundary thus acts as a barrier for the ions. Consequently, a null flux is imposed 
(natural boundary conditions). 

[0142] The imposition of the boundary conditions for the water content 0 is more complicated. It 
was mentioned that the soil S2 is at an 85% relative humidity and that the environment 
50 is at 70%. To convert these values in water content, the water adsorption isotherm 
model is used, which gives the nonlinear relationship between 0 and the relative 
humidity in the material (see Figure 6). According to this model, an 85% relative 
n humidity corresponds to a 0.0823 m 3 /m 3 water content, and a 70% relative humidity 
■0 corresponds to 0.0626 m 3 /m 3 . These values are imposed at x=0 and x=0.15 
ill respectively. 

[0fi§3] For the electrical potential % a reference value must be set at some point. A value of 
zero is thus imposed at x=0. 

[Ofy.4] The initial concentration in the pore solution of the material is given in Table 2. It was 
Q evaluated with the chemical equilibrium code developed by Reardon in "An ion 
5=58 interaction model for the determination of chemical equilibria in cement/water systems", 
Cement and Concrete Research, vol. 20, p. 175-192, 1990. The initial value of the water 
content is 0.0858 m 3 /m 3 , which corresponds to a relative humidity of 87% in the 
material. Finally, the initial value of the potential is zero throughout the slab. 

[0145] Table 2 : Data for the numerical simulations - First example 



Properties 


Value 


Properties 


Value 


Cement type 


V 


Cement composition 


(%) 


w/c 


0.50 


C 3 S 


56.31 






C 2 S 


22.04 






C3A 


1.69 
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C4AF 


11.26 


Concrete mixture 


(kg/m J ) 


Initial solid phases 


(g/kg) 


Cement 


342.5 


C-S-H 


89.4 


Water 


171.3 


Portlandite 


44.8 


Aggregates 


1944.0 


Ettringite 


14.6 


Density 


2457.8 


Hydrogarnet 


11.1 


Initial pore solution 


(mmol/L) 


Diffusion coefficients 


(m z /s) 


OH" 


626.0 


OH" 


2.03x10" 10 


Na + 


204.1 


Na + 


5.13x10' 11 


K + 


475.8 


K 


7.53x10' 11 


S0 4 2 ' 


28.0 


S0 4 2 " 


4.10x10' 11 


Ca 2+ 


1.1 


Ca 2+ 


3.05x1 0' 11 


AI(OH) 4 " 


0.1 


AI(OH) 4 " 


2.08x10" 11 


Rel. humidity (%) 


(rrv7m°) 


Water diffusivity 


(m z /s) 


70 


0.0626 




2.72x1 0" 12 e 81 - 29 


85 


0.0823 






87 


0.0858 






Porosity 


0.118 


R 


0.0 \ho j/moi/ l\ 


T 


25°C 


£ 


6.94x1 0" lu CA//m 


F 


96488.46 C/mol 


T 


0.038 


e ° 


1.602x10"* C 







[0146] All the other parameters and characteristics of the material are given in Table 2. The 
initial distribution of solid phases is calculated by considering the hydration reactions of 
cement. It is assumed that the paste is fully hydrated at the beginning of the 
simulations. The porosity is evaluated from data on hydrated cement pastes found in 
Taylor H.F.W., "Cement chemistry", Academic Press (Canada), 1990. The diffusion 
coefficient values are based on results of chloride migration experiments. These values 
are used to evaluate the tortuosity t, which is given by the diffusion coefficient of a given 
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species in the material (D/) divided by its value in free water {DC). Taking for example 
OH", D, is 2.03x10' 10 m 2 /s and DC is 5.273x10" 9 m 2 /s, which yields a tortuosity of 0.038. 
The adsorption of ions is neglected in this example. 

[0147] The only parameters not yet discussed related to the movement of water in the pore 
system under capillary effects. As seen in equations (7) and (14), two parameters are 
needed: D e and D L . However, it was shown that for high values of relative humidity, 
such as in the present example, D e * D L . Hence, only the value of D L is needed. Based 
on NMR imaging of water absorption in mortars, this nonlinear parameter is taken as 
q D L =2.72x10" 12 e 81 - 29 . 

[0118] Before showing any simulation results, note on the ability of the chosen algorithm to 
S handle the water flowing through the material as a result of capillary suction. To have an 
5 idea of the relative importance of the flow with regards to diffusion, the Peclet number 
(P e ) for the material properties mentioned previously is calculated. P e is given by VL/D, 
U where V is the velocity of the flow, L is the length of the domain and D is the diffusion 
^ coefficient. In the present case, V can be calculated with equation (8), with the gradient 
5 of water content evaluated from the analytical solution of equation (7) in steady-state. It 
|aS gives 1.42x10" 10 m/s, corresponding to a water flux of 12.3 ml/day/m 2 . Since multiple 
ions are considered, one species must be selected in order to calculate P e . The 
diffusion coefficient of OH", which is 2.03x1 0" 10 m 2 /s is chosen. This results in a Peclet 
number of 0.11, which is low. The problems commonly associated with convection- 
dominated problems, like strong oscillations, are not likely to occur at such low values of 
P e . Even for a very high gradient of relative humidity such as 95% - 60%, P e is 0.7, still 
well below values expected to lead to numerical oscillations. 

[0149] But still, two types of problem occur. For the particular case of a concrete slab S2 
exposed on its upper surface to the environment, the natural boundary condition at this 
location prevents the ions from exiting the slab SO. Combined with the upward water 
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flux, it means that the ions tend to accumulate in the upper part of the slab. In the case 
of a sulfate attack, this is particularly true for Na + and S0 4 2 " since they are present in the 
groundwater. This can lead to very high concentration gradients at x=L, and 
consequently to divergence of the algorithm. Numerical tests have shown that at some 
point, the concentrations in Na + and S0 4 2 " may reach such high levels that even 
increasing the number of elements in this region cannot prevent the algorithm to crash. 
Taking the formation of mirabilite into account in the model helps to prevent this 
situation since it absorbs a large amount of Na + and S0 4 2 " ions. 

[0150] The other problem caused by the presence of moisture transport in the model occurs in 
2 the first time steps. Starting from an initial water content level, imposing Dirichlet-type 
1 conditions may lead to numerical oscillations near the boundaries because of initially 
w high gradients in these areas. In a classical diffusion problem, these oscillations would 
I tend to disappear with time. With the current model, however, they must be avoided. 
F Since the concentration is coupled to the water content, oscillations in the latter can 
□ induce oscillations in the ionic profiles. The problem is that if the oscillations are such 
nl that they lead to negative concentrations at some nodes, the chemical equilibrium 
£ module will not be able to find a solution since negative concentration do not make 
P physical sense. To avoid this problem, one solution consist in increasing the number of 
nodes near the boundaries, to avoid, as much as possible, the oscillations. Another 
solution is to gradually impose the boundary conditions for the water content over time, 
which is similar to the incremental load method commonly used in solid mechanics. This 
technique proved successful to avoid the problems related to the oscillations in the 
water content profiles. Finally, if the initial humidity in the material is lower than 100%, 
oscillations may not be a problem. With the data proposed in Table 2, i.e. a 87% RH 
initial material exposed to a 85% - 70% RH gradient, none of the previously exposed 
problems happens. 

[0151] In this example, the simulations are performed over a period of twenty years. The 
domain is discretized with 50 elements of equal length. The time steps are one day 
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long. 

Results and discussion 

[0152] The solid phase distribution is shown in Figure 7. One of the main effect of the 
penetration of sulfates in the material is the formation of an ettringite front that proceeds 
toward the top of the slab. It is located at 3 cm from the exposed surface after 20 years. 
It is accompanied by a narrow peak of gypsum that forms near the bottom part of the 
slab, between 0 and 1 cm. Since calcium ions are necessary to form ettringite, it can be 
seen that a small portion of portlandite, whose position matches with the front of 
*0 ettringite, was dissolved in order to achieve that. Also matching the position of the 
in ettringite front is the dissolution of hydrogarnet, whose Ca 2+ and AI(OH) 4 ~ are used to 
^ form ettringite. When the hydrogarnet is totally dissolved, no further amount of alumina 
O is available, which explains that the formation of ettringite reaches a maximum. It must 
be noted that ettringite also precipitates near the upper surface of the slab. Finally, 
y portlandite is completely dissolved between the exposed surface and 0.5 cm, while only 
fU a small amount of C-S-H is decalcified. 

[0V53] Since it is assumed that there are no Ca 2+ and OH" ions in the soil, they tend to leach 
out of the material. This leads to a dissolution of portlandite and a decalcification of C-S- 
H. Again, these chemical reactions tend to progress upward. The numerical results 
clearly show that the dissolution of portlandite occurs before the decalcification of C-S- 
H, an observation also made in experimental tests. 

[0154] Finally the results show the formation of mirabilite on the upper part of the material. It is 
caused by the accumulation of Na + and S0 4 2 " ions in this region, as a result of the flow 
of water directed upward. This result makes sense since the upper part of the slab is the 
only place in the material where the Na + and SO4 2 " ions can accumulate and reach the 
concentration level necessary to precipitate. The precipitation of ettringite at this 
location may also be an effect of the upward water flow. 
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[0155] Figure 8 shows the distribution of ions after ten years. The profile of AI(OH) 4 " is not 
shown since its concentration values are very small. A verification made after the 
calculation showed that the ionic concentrations at a given nodes are in equilibrium with 
the solid phases found at this particular location. 

[0156] Simulations were performed without the electrical coupling and the chemical activity 
effect to weigh their effect on the overall result. The comparison with the previous result 
is shown on Figure 9. Neglecting the chemical activity effects does not seem to have an 
important effect on the ettringite front. Yet, the model does not predict the formation of 
5 ettringite at the surface (x=0.15) when the influence of chemical activity is not 
1 considered. So obviously, chemical activity is an important part of the process and 
4? should be considered. When the electrical coupling is also neglected, the difference in 
3 the position of the ettringite front is important. This shows the importance of considering 
* both the electrical coupling and the chemical activity. 

[0ji7] The next results show the effect of the capillary suction on the solid phase distribution. 
£ Figure 10 compares the ettringite and portlandite profiles for the saturated and 
P unsaturated cases. The simulation for the saturated case was performed with an initial 
water content in the material that is equal to the porosity, meaning that all of the porous 
space is filled with water. The boundary conditions for the water content also 
corresponds to the porosity for both ends of the slab. All of the other boundary 
conditions as well as all of the other parameters are the same as for the previous 
simulations. The results show that the ettringite fronts penetrates less in the material 
when it is saturated. It clearly emphasizes the effect of the upward water flow caused by 
capillary suction on the sulfate ions, that go farther in the concrete in that case. The 
difference between the portlandite profiles is very slight. Nevertheless, one can see that 
the degradation of portlandite is less important for the unsaturated case. It shows that 
the upward flow of water tends to slow the downward leaching of Ca 2+ and OH" ions, 
thus slowing the dissolution of portlandite. Simulations have shown that the effects on 
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both the ettringite and the portlandite described above are more important as time 
increases. 

[0158] It was mentioned earlier that the algorithm needs small time steps in order to perform 
properly. The next simulations were performed to measure the effect of different time 
steps on the solution. Figure 11 shows the ettringite profile after twenty years obtained 
through time steps of 0.5, 1, 2, 5, and 10 days. It shows that as the time steps is 
increased, the ettringite front is less and less sharp. Thus, the time step clearly has an 
influence on the shape of the profile, with a smaller time step meaning sharper, more 
precise solid phase profiles. As the time steps are getting smaller, the difference 
S between each solution reduces, ultimately converging toward a unique solution. In 
£ Figure 11, the difference between the solution with At="\ day and At=0.5 day is very 
W small, so the one day time step is sufficient to reliably perform the simulation. 

[0j$9] The same kind of simulation was performed but this time with a different number of 
O elements. The results are given in Figure 12 for 30, 50, 100 and 200 elements, again for 
ill the ettringite phase. The profiles show some differences near the boundary x=0. 
J However, this is simply a discretization effect. For all cases, there is no ettringite on the 
S=^ first node, and it reaches its maximum value on the next node. What is more interesting 
is the solution at the position of the ettringite front. The solution is almost identical, 
whatever the number of elements, except maybe for the case with 30 elements. Even in 
that case, the difference is not very important compared to those obtained with different 
time steps. It means that the choice of element should be made in order for the 
transportation part of the algorithm to converge. This choice has virtually no influence 
on the shape of the solid phase profiles. 

SECOND EXAMPLE 

Description of the problem 
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[0160] The model can be used to predict the degradation of various cement-based materials 
(pastes, mortars, concretes) exposed to aggressive chemical solutions. It is used herein 
to analyze the behavior of cement paste disks made of an ordinary Portland cement 
prepared at a water/cement ratio of 0.60. The paste was cured for a year in aluminum 
foil. It was then cut in 12 mm disks. As shown in Figure 13 these samples 57 were 
immersed for three months in a sodium sulfate solution 55 prepared at a concentration 
of 50 mmol/L. All sides of the sample except one were coated with silicon 58. The 
material being exposed to the aggressive solution on only one face ensures a one- 
dimensional process. In order to maintain constant the boundary conditions during the 
test, the test solution was regularly renewed. The microstructural alterations were 
3 investigated by means of electron microprobe analyses and scanning electron 
microscopy. 

[0O1] The various characteristics of the material are given in Table 3. The initial solid phases 
=p forming the cement paste are calculated by considering the hydration reactions of 
O cement. The results are adjusted to take into account the one-year curing time. All of 
m the S0 3 present in the cement reacts with the C3A to form the initial amount of ettringite. 
:£ The remaining part of the C3A as well as the C4AF react with water to yield hydrogarnet. 
H Following microscope observations, it is assumed that only 35% of C4AF reacts. 

[0162] The initial pore solution is obtained according to the pore solution extraction procedure. 
Samples were placed in an extraction cell and crushed at a pressure of approximately 
300 MPa. Typically, 2 to 5 ml of pore solution are extracted. The solution is delivered 
through a drain ring and channel, and recovered with a syringe in order to limit exposure 
to the atmosphere. Chemical analyses of the pore solution were carried out shortly after 
the extraction test. The composition of the solution was then adjusted in order to respect 
the electroneutrality condition. Slight deviations from electroneutrality can arise due to 
the experimental error associated with the evaluation of the ionic concentrations. The 
neutral solution is once again corrected so that the concentrations of the species are in 
accordance with the chemical equilibrium constants displayed in Table 1. This 
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correction is performed with the chemical equilibrium code described previously. 

[0163] The porosity is calculated according to the standardized ASTM C642 water porosity 
procedure. 

[0164] The diffusion coefficient of each ionic species is evaluated according to the procedure 
detailed in the second embodiment of the invention. The upstream and the downstream 
cells of the migration setup both contain a 300 mmol/L NaOH solution. The upstream 
cell also contains a 500 mmol/L NaCI solution. Two disks are tested. They are about 25 
mm thick with an exposed surface of approximately 75 mm. They were both vacuum 
S saturated for 24 hours prior to the test. A 5 Volt potential difference is applied to the 
!| sample for the full duration of the test. Electrical current through the setup is measured 
W for four days (90 hours). A preferred embodiment of the invention is then used to 
3 reproduce the measured currents. The thickness of the sample is discretized with 200 
^ elements. The mesh is refined in a 2 mm layer on both extremities of the domain, where 
O 50 elements are used. The simulations with various values of t are performed with 90 
m time steps of one hour each. The value of x that offers the best match with the 
f measured current curve yields the diffusion coefficients for that particular material. This 
best-fitting curve is shown in Figure 14 for one of the disks. A similar curve is obtained 
for the other disk. The diffusion coefficients shown in Table 3 correspond to the average 
values for the two sets of measurements. The corresponding tortuosity value is also 
indicated in Table 3. 

[0165] This example occurs in saturated conditions. Consequently, no capillary suction occurs 
in the materials. This is modeled by setting D&=D L =0. 

[0166] All of the parameters needed to start the simulations are now determined. They will 
serve as input parameters to the preferred embodiment of the invention. The 
simulations are performed over a three-month period using 1080 time steps of two 
hours. The 12 mm disks are discretized with 50 elements. The latter all have the same 
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length. 



[0167] Table 3 : Data for the numerical simulations - Second example 
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Hvdroaarnet 


39.2 


Density 


1744.0 






Initial pore solution 


(mmol/L) 


Diffusion coefficients 


(m z /s) 


OhT 


433.2 


OH" 


1.61x10' 10 


Na + 


111.1 


Na + 


4.07x10" 11 


K + 


327.0 


K + 


5.98x1 0' 11 


S0 4 2 ~ 


4.0 


S0 4 2 ~ 


3.25x10- 11 


Ca 2+ 


1.6 


Ca 2+ 


2.42x10" 11 


AI(OH) 4 ' 


0.07 


AI(0H) 4 " 


1.62x10' 11 


Porosity 


0.522 


R 


8.3143 J/mol/°K 


T 


25°C 


s 


6.94x10- 1u CA//m 


F 


96488.46 C/mol 


X 


0.038 


e 0 


1.602x10" ia C 







Results and discussion 
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[0168] The invention predicts the spatial distribution of the various solid phases after a given 
period of exposure. The results for the 0.60 paste mixture immersed for 3 months in the 
sodium sulfate solution are summarized in Figure 15. This ability of the model to predict 
the distribution of various solid phases should be of particular interest for those 
concerned by the consequences of chemical damage on the mechanical properties of 
the material. The precise knowledge of the solid phase distribution allows to calculate 
(at each location) the actual porosity of the material. These results could eventually be 
used to estimate the residual mechanical strength of the material. 

[01J>9] The microprobe analysis gives profiles of total calcium and total sulphur content. To 

3 make comparisons with the results given by the invention, the total calcium content is 

J2 calculated by considering the calcium found in portlandite, C-S-H, ettringite, hydrogarnet 

W and gypsum. This gives the profile of Figure 16. The same procedure is applied to 

3 sulphur, which is found in ettringite and gypsum. The result is shown in figure 17. 

[OPO] As can be seen in Figure 16, the algorithm of the preferred embodiment matches very 
ill well with the measured calcium profile. The invention can predict very accurately the 
f beginning of the calcium loss, located at 1.3 mm from the surface. Toward the core of 
P the sample (increasing x), no sign of degradation is apparent, which the model also 
reproduces. This is a very important result since the loss of calcium comes from the 
dissolution of portlandite and the decalcification of C-S-H, the two solid phases that give 
the material its strength. The invention can thus be used to predict the resistance loss of 
cementitious materials undergoing calcium leaching. 

[0171] Good results are also obtained with regard to the sulfur profile (Figure 17). The 
microprobe analysis shows a slight increase in sulfur above the base level in the core of 
the sample at 2.2 mm. Not only does the model matches this increase but it gives 
information on the nature of this phenomenon. By comparing Figure 17 with Figure 15, 
one can see that it can be attributed to the penetration of an ettringite front. 
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[0172] The invention predicts a peak of sulfur near the surface of the material, also observed 
with the microprobe analysis. The comparison with Figure 15 reveals that it is a peak of 
gypsum. There is a slight offset between the location of the peak measured with the 
microprobe and the one predicted by the model. The important amount of gypsum 
precipitated near the surface causes a very important reduction of the porosity in this 
region. Accordingly, the invention, through the damage factor G (see equation 20), 
reduces the diffusion coefficient of each ionic species. However, in the paste, the 
amount of gypsum is so important that it causes microcracking. This increases the 
transport properties instead of lowering them. It explains that the peak predicted by the 
model lags behind the one shown by the microprobe analysis. 

It should be noted that the present invention can be carried out as a method, can be 
embodied in a system, a computer readable medium or an electrical or electro- 
magnetical signal. 

[0f?4] It will be understood that numerous modifications thereto will appear to those skilled in 
"Si the art. Accordingly, the above description and accompanying drawings should be taken 
£ as illustrative of the invention and not in a limiting sense. It will further be understood 
P that it is intended to cover any variations, uses, or adaptations of the invention following, 
in general, the principles of the invention and including such departures from the 
present disclosure as come within known or customary practice within the art to which 
the invention pertains and as may be applied to the essential features herein before set 
forth, and as follows in the scope of the appended claims. 
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